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Eilenberger equations, and use a variation of the Brandt-Pesch-Tewordt (BPT) method to obtain a 
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two-dimensional dj.2_y2 {dxy) superconductors with the field rotated in the basal plane. The solution 
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the gap nodes may result in minima or maxima in the angle-dependent specific heat, depending on 
the location in the T-H plane. This variation is attributed to the scattering of the quasiparticles on 
vortices, which depends on both the field and the quasiparticle energy, and is beyond the reach of 
the semiclassical approximation. We investigate the anisotropy across the T-H phase diagram, and 
compare our results with the experiments on heavy fermion CeCoIus. 
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I. INTRODUCTION 

In this paper and its companion^, hereafter referred to 
as II, we present a general theoretical approach for in- 
vestigation of thermal and transport properties of super- 
conductors in magnetic field, and use it to determine the 
behavior of the density of states, specific heat, and ther- 
mal conductivity in the vortex state of unconventional 
superconductors. Our more specific goal here is to pro- 
vide connection between theory and recent experiments 
measuring the properties of such superconductors under 
a rotating magnetic field, to explain the existing data, 
and to guide future experimental studies. We focus on 
these experiments as they hold exceptional promise for 
helping determine the structure of the superconducting 
energy gap. 

We consider unconventional superconductors, for 
which in the ordered state both the gauge symmetry and 
the spatial point group symmetry are broken^. Then the 
gap in the single particle spectrum, |A(p)|, is momentum 
dependent. We focus on anisotropic pairing states with 
zeroes, or nodes, of the superconducting gap for some 
directions on the Fermi surface (FS). 

The single particle energy spectrum of a superconduc- 
tor is E{ii) = ±^^2(p) ^ |A(p)|2, where ^(p) is the 
band energy in the normal state with respect to the Fermi 
level. Consequently, the gap nodes, |A(p)| = 0, arc the 
loci of the low energy quasiparticles, and the number of 
quasiparticles excited by temperature or other pertur- 
bations depends on the topology of the nodal regions. 
Experimental probes that predominantly couple to un- 
paired electrons, for example the heat capacity or (for 
pairing in the singlet channel) magnetization, are com- 
monly used to show the existence of the gap nodes. The 
nodal behavior is manifested by T" power laws, with the 
exponent n that depends on the structure of the gap^. 

Locating the nodes on the Fermi surface is a harder 



task. Since usually only the phase of the gap, but not the 
gap amplitude, |A(p)|, breaks the point group symmetry, 
transport coefficients in the superconducting state retain 
the symmetry of the normal metal above Tc- The phase 
of the order parameter can be tested by surface mea- 
surements, but experimental determination of the nodal 
directions in the bulk requires breaking of an additional 
symmetry. One possible approach is to apply a magnetic 
field, H, and rotate it with respect to the crystal lattice. 
The effect of H on the nodal quasiparticles depends on 
the angle between the Fermi velocity at the nodes and 
the field, and hence provides a directional probe of the 
nodal properties^. 

At the simplest level, screening of the field and the re- 
sulting flow of the Cooper pairs, either in the Meissner or 
in the vortex state, locally shifts the energy required to 
create an unpaired quasiparticle relative to the conden- 
sate (Dopplcr shift )iAs§.. Our focus here is on the vortex 
state, where the supercurrents are in the plane normal to 
the applied field, and hence only the quasiparticles mov- 
ing in the same plane are significantly affected. Applying 
the field at different angles with respect to the nodes pref- 
erentially excites quasiparticles at different locations at 
the Fermi surface, and leads to features in the density 
of states (as a function of the field direction)^. This, in 
turn, produces oscillations in the measurable thermody- 
namic and transport quantities, which can be used to 
investigate the nodal structure of unconventional super- 
conductors. 

Such investigations have been carried out exper- 
imentally in a wide variety of systems. Due to 
higher precision of transport measurements, more 
data exist on the thermal conductivity anisotropy 
under rotated field. The anisotropy was reported 
in high-temperature superconductors^!^, heavy fermion 
UPd2Al3S, CcCoIngia, PrOs4Sbi2^, organic k-(BEDT- 
TTF)2Cu(NCS)2^, and borocarbide (Y,Lu)Ni2B2GH, 
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see Ref. [13 for review. The heat capacity measure- 
ments are more chahenging, and were carried out in 
the borocarbidesi^^, and CeCoIng^. While the experi- 
ments provided strong indications for particular symme- 
tries of the superconducting gap in these materials, they 
did not lead to a general consensus. The main reason 
for that has been lack of reliable theoretical analysis of 
thermal and transport properties in the vortex state. 

Historically, there was a schism between theoretical 
studies of the properties of s-wave type-II superconduc- 
tors at low fields, where the single particle states are 
localized in the vortex cores, and the investigations near 
the upper critical field, i?c2, where vortices nearly over- 
lap and the quasiparticles exist everywhere in space. The 
distinction between the two regimes is not so clear cut in 
unconventional superconductors, since it is the extended 
near-nodal states that control the electronic properties 
both at high and at low fields^. Often it is hoped that 
a single theoretical approach may provide results valid 
over a wide temperature and field range in nodal super- 
conductors. 

In part because early experiments on the vortex state 
of unconventional superconductors focused on the high- 
Tc cuprates^i^, theoretical work has long been rooted in 
the low field analysis. The Dopplcr shift approximation 
was used to predict and analyze the behavior of the spe- 
cific heat^i^ii^ and the thermal conductivityi^iSSiSii^ un- 
der an applied magnetic field. The method is semiclassi- 
cal in that it considers the energy shift of the nodal quasi- 
particles with momentum p at a point R. Consequently, 
it is only valid at low fields, H <C Hc2, when the vortices 
are far apart, and the supervelocity varies slowly on the 
scale of the coherence length. Moreover, most such cal- 
culations account only for quasiparticles near the nodes, 
and therefore are restricted to energies small compared 
to the maximal superconducting gap, and hence to tem- 
peratures T <C Tc- In addition, the energy shift leaves 
the quasiparticlc lifetime infinite in the absence of impu- 
rities, and therefore the method docs not account for the 
scattering of the electrons on vortices. While some at- 
tempts to remedy the situation existii^i^, no consistent 
description emerged. 

Recent experiments cover heavy fermion and other low 
temperature superconductors, and generally include the 
regime T < and H < Hc2- Consequently, there has 
been significant interest in developing alternatives to the 
low field Doppler shift approach. The goal is to treat 
transport and thermodynamics on equal footing, to be 
able to describe the electronic properties over a wide 
range of fields and temperatures, and to include the ef- 
fects of scattering on vortices. Fully numerical solution 
of the microscopic Bogoliubov-dc Genncs equations have 
been employed for computing the density of states (see, 
for example, Ref. [23 ). but are not naturally suited for 
computing correlation functions and transport proper- 
ties. Calculation of the Green's function in the super- 
conducting vortex state is difficult due to appearance of 
additional phase factors from the applied field. More- 



over, transport calculations need to include the vertex 
corrections, since the characteristic intervortex distance 
is large compared to lattice spacing, hence the scattering 
on the vortices corresponds to small momentum transfer, 
and the forward scattering is important. 

Here we use the microscopic approach in conjunction 
with a variant of the approximation originally due to 
Brandt, Pesch and Tewordt (BPT)^ that replaces the 
normal electron part of the matrix Green's function by 
its spatial average over a unit cell of the vortex lattice. 
While originally developed for s-wave superconductors, 
this approach has recently been successfully and widely 
applied to unconventional systems (see Sections IIIBI and 
mil for full discussion and references), where it gave re- 
sults that are believed to be valid over a wide range of 
temperatures and fields^^iSl. 

We employ the approximation in the framework of the 
quasiclassical metho d^^i^^ . Two main advantages of this 
approach are: a) BPT approximation results in a closed- 
form solution for the Green's function2£i^i^ enabling us 
to enforce self-consistency for any field, temperature, and 
impurity scattering, and facilitating the subsequent cal- 
culations of physical properties; b) quasiclassical equa- 
tions are transport-like, so that the difference between 
single particle and transport lifetimes appears naturally, 
without the need to evaluate vertex corrections. Con- 
sequently, we are able to compute the density of states, 
specific heat, and the thermal conductivity on equal foot- 
ing, and provide a detailed comparison with experiment. 

In this we pay particular attention to the data on heavy 
fermion CeCoIns, where the specific heat and the ther- 
mal conductivity data were interpreted as giving contra- 
dictory results for the shape of the superconducting gap. 
The anisotropic contribution to the specific heat exhib- 
ited minima for the field along the [100] directions, which 
led the authors to infer d^y gap symmetryii, while the 
(more complicated) pattern in the thermal conductivity 
for the heat current along the [100] direction under ro- 
tated field was interpreted as consistent with the d^^-y^ 
gap>iS. In a recent Letter— we suggested a resolution for 
the discrepancy, and provide the detailed analysis here. 

The remainder of the paper is organized as follows. 
In Sec. [ni we briefly review the quasiclassical approach 
and the BPT approximation to the vortex state. Sec. IIIII 
gives the derivation of the equilibrium Green's function. 
Some of the more technical aspects of the calculation are 
described in the appendices: Appendix \X\ describes a use- 
ful choice of ladder operators that enable us to efficiently 
solve the quasiclassical equations in the BPT approach, 
and Appendix [Blshows how to find a closed form solution 
for the Green's function. 

Many of the salient features of our results are clear 
from a simple and pedagogical example of a 2D rf-wave 
superconductor with a cylindrical Fermi surface consid- 
ered in Sec. IIVBI We discuss the influence of the field 
on the density of states in the vortex state and present 
the results for the anisotropy of the specific heat and 
heat conductivity for an arbitrary direction of the ap- 
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FIG. 1: (Color online) In this paper we present calculations 
for a simple cylindrical Fermi surface and a Fermi surface 
shown here. The d-wave order parameter has lines of vertical 
nodes. Our goal is a calculation of the thermodynamic prop- 
erties, such as specific heat and entropy, and their anisotropy 
under magnetic field rotations, (po, in the a&-plane. 

plied magnetic field. 

As one of our goals is the comparison of the results 
with the data on layered CeCoIns, Sec. IIV CI is devoted 
to the fully self-consistent calculations for more realistic 
quasi-cylindrical Fermi surfaces, Fig.[T] The discussion of 
the results, comparison with the data, and implications 
for future experiments are contained in Scc.|Vl The com- 
panion paper II uses these results to derive and discuss 
the behavior of the thermal conductivity. 

We aimed to make the article useful to both theorists 
and experimentalists. Sec. IIVBI and Sec. |V] arc proba- 
bly most useful for those readers who are interested only 
in the overall physical picture and the behavior of the 
measured properties; the figures in Sec. IIV CI show the 
main differences between the self-consistent and non-self- 
consistent calculations. 



II. QUASICLASSICAL APPROACH 

A. Basic equations and formulation 

We begin by writing down the quasiclassical 
equations for a singlet superconductor in magnetic 
gg y28, 29^31, 33^34, 35 summarizing the details relevant 
for our discussion. The equations are for the quasiclassi- 
cal (low-energy e) Green's function, which is a matrix in 
the Nambu (spin and particlc-hole) space, 

J(R,p;.,.(4-/)^ ,1, 

This matrix propagator has been integrated over the 
quasiparticle band energy, and therefore depends only 
on the direction at the Fermi surface, p, and the center 
of mass coordinate, R. 



We formulate our approach in terms of the real-energy 
retarded, advanced, and Keldysh propagators. This is 
a natural path for the self-consistent calculation of the 
quasiparticle spectrum, needed for determination of ther- 
modynamic properties such as entropy and heat capac- 
ity. Moreover, the Keldysh technique is the most direct 
route towards non-equilibrium calculations, required for 
the transport properties such as thermal conductivity, 
which is covered in the companion paper II. Consequently 
we establish a unified approach to describe both the ther- 
modynamics and transport in the vortex state. 

Retarded (R) and advanced (A) functions g = g^'^ 
satisfy (we take the electron charge e < 0) 

[(£ + ^v/(p)A(R)) ?3 - A(R, p) - ct™p(R; e), 

5(R, p; e)] + iv/(p) • V^, g{R, p; e) = , (2) 
together with the normalization condition 

?^-^(R,p;£„,)2 = _^2^. (3) 

Here e is the real frequency, v/(p) is the Fermi velocity 
at a point p on the FS. The magnetic field is described by 
the vector potential A(R), and the self-energy a is due 
to impurity scattering. The equations for the retarded 
and the advanced functions differ in the definition of the 
corresponding self-energies. 

The mean field order parameter, 

is defined via the self-consistency equation involving the 
Keldysh function f^, 

A(R,p)= rfp^s n/(p')V(p,p')/'^(R,P';£), 

(5) 

In equilibrium /^^ = (/^ — /'^) tanh(e/2T), and we ob- 
tain the usual self-consistency equation computing the 
£-integral in the upper (lower) half-plane for (/^)- 

We wrote Eq.fO for a general Fermi surface, and there- 
fore introduced the density of states (DOS) at a point p 
on the Fermi surface in the normal state, Nf{p). The 
net density of states, Nf = / dppsNf{p), and we define 
n/(p) = Nf{f))/J\ff. We absorbed the net DOS, TV/, into 
the definition of the pairing potential, V{p,p'). 

Since below we frequently perform the integrals over 
the Fermi surface, we introduce a shorthand notation 

( • )fS = J rfPFS «/(p) • , (6) 

SO that the the gap equation above can be rewritten as 
A(R, p) = 1 ^ ( V{p, p') /^^(R, p'; e))ps ■ (7) 
All calculations below arc for separable pairing, 

v{p,p')^v,y{p)y{p'), (8) 
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where 3^(p) is the normahzcd basis function for the par- 
ticular angular momentum, (3^(p)^)fs = 1- For example, 
for da;2_y2 gap over a Fermi surface parameterized by 
angle we have y{(j)) = \/2cos20. Hence the order 
parameter is A(R, p) = A(R)3^(p). 

Finally, we include the isotropic impurity scattering 
via the self-energy. 



(9) 

Here riimp is the impurity concentration, and, in the self- 
consistent ^-matrix approximation, 

<(R;e) = ui + u7V/(g(R,p;e))ps <(R;e), (10) 



where u is the single impurity isotropic potential. Com- 
paring Eq. ([2]) and Eq. ([9]) we see that S effectively rcnor- 
malizes the energy e, while A^mp accounts for the impu- 
rity scattering in the off-diagonal channel. The term D 1 
drops out of equations for the retarded and advanced 
Green's functions since the unit matrix commutes with 
the Green's function in Eq. ([2|). This term, however, 
generally appears in the Kcldysh part, and has a sub- 
stantial effect on transport propertieSi ^^'^^ Below we pa- 
rameterize the scattering by the "bare" scattering rate, 
r = nimp/i^Nf, and the phase shift i5o of the impurity 
scattering, tani^o = ttuA//. 



In equilibrium wc explicitly write Eqs. (H))-® as a system of equations, 

»v/(p)- V^.g + A/-A/ = 0, 

/ Oip \ 1 ~ 

/ = 2*Ag, 
./ = 2zA5, 



-2z£-fv;(p) ( v^-^a(r; 



-2ie-v;(p) V 



2ie 



A(R) 



(11a) 
(lib) 

(11c) 
(lid) 



where e = e— E, A = A 



Ajmp, and A = A* 



B. Vortex state ansatz and Brandt- Pesch-Tewordt 
approximation 

So far our discussion remained completely general. In 
the vortex state of a superconductor, the order parameter 
and the field vary in space, and the quasiclassical equa- 
tions have to be solved together with the self-consistency 
equations for the gap function, and Maxwell's equation 
for the self-consistently determined magnetic field and 
the vector potential. Finding a general non-uniform so- 
lution of such a system is a daunting, or even altogether 
impossible, task. Therefore we make several simplifying 
assumptions and approximations that allow us to obtain 
a closed form solution for the Green's function. 

First, wc assume the magnetic field to be uniform. This 
assumption is valid for fields H ^ Hd , where the typical 
intervortex spacing (of the order of the magnetic length, 
A = (fic/2|e|i3)^/^) is much smaller than the penetration 
depth, the diamagnetic magnetization due to the vor- 
tices is negligible compared to the applied field, and the 
local field is close to the applied external field, B k, H . 
All the materials for which the anisotropy measurements 
have been performed are extreme type-II superconduc- 
tors, where this assumption is valid over essentially the 
entire field range below Hc2- 

In writing the quasiclassical equations we only included 
the orbital coupling to the magnetic field, assuming that 



it dominates over the the paramagnetic (Zccman) con- 
tribution. This is valid for most superconductors of in- 
terest, and the detailed analysis of the Zeeman splitting 
will be presented separately — ;the main conclusions of 
this paper remain unaffected. 

Second, we take an Abrikosov-like vortex lattice ansatz 
for the spatial variation of the order parameter, which is 
a linear superposition of the single- vortex solutions in the 
plane normal to the field. We enforce the self-consistency 
condition, which requires going beyond the simple form 
suggested by the linearized Ginzburg-Landau equations. 
The details of this choice arc given in Sec lIIII below. 

In the vortex lattice state the quasiclassical equations 
generally do not allow solution in a closed form. We 
therefore employ a variant of the approximation origi- 
nally due to Brandt, Pesch and Tewordt (BPT)^. The 
method consists of replacing the diagonal part of the 
Green's function by its spatial average, while keeping 
the full spatial structure of the off-diagonal terms. It 
was initially developed to describe superconductors near 
the upper critical field, where the amplitude of the or- 
der parameter is suppressed throughout the bulk, and 
the approximation is nearly exact. This is confirmed 
by expanding the Green's function in the the Fourier 
components of the reciprocal vortex lattice, (?(R, p; s) = 
(7(K, p; e) exp(iKR). and noticing that g^(K.) oc 
exp(— A^K^) so that the K = component is expo- 
nentially dominant In situations where the states in- 
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side vortex cores are not crucial for the analysis, such 
as in extreme type-II or nodal superconductors^Si^ 
the method remains valid essentially over the entire field 
range. Consequently the BPT approach and its varia- 
tions was extensively used to study unconventional su- 
perconductors in the vortex statej^i^>^ One of the ad- 
vantages of the method that it reproduces correctly the 
H = BCS limit)^ and therefore may be used to inter- 
polate over all fields. One, however, needs to be cautious 
in computing the properties of impure systems: averag- 
ing over the intervortex distance A) prior to averag- 
ing over impurities is allowed only when A/£ ^ 1, where 
£ is the mean free path, and hence the approach does 
break down at very low fields, and only asymptotically 
approaches the zero field result. We show the signatures 
of this breakdown in Sec. IIV CI 

The use of the BPT approximation relaxes the con- 
straints imposed by the assumption of a perfectly pe- 
riodic vortex arrangement. Indeed, averaging over the 
unit cell of the vortex lattice is somewhat akin to the 
coherent potential approximation in many body physics, 
although with an important caveat that this is only done 
for the normal part of the matrix Green's function. Con- 
sequently, the results derived within this approach arc 
also applicable to moderately disordered vortex solids. 

III. SINGLE-PARTICLE GREEN'S FUNCTION 

Hereafter we use g to denote the spatially av- 
eraged electron Green's function, g = g(p; e) = 
(7(R, p;e). The approach we take here follows the stan- 
dard practice^i^i^'^ of determining g from the spatially 
averaged normalization condition, Eq. (|llap . 

5' -77 = -^'- (12) 

Here we defined the average over vortex lattice of a prod- 
uct as 

/HTt 
— /i(R)/2(R). (13) 

The anomalous components of the Green's function sat- 
isfy Eqs. (|llc]l - (jlldp . Formally, the solution is obtained 
by acting with the inverse of the differential operator in 
the right hand side on the product A g and A g respec- 
tively. Upon replacement of g by its average, the operator 
acts solely on the order parameter, 

/(R,p;e)= 2i5(p;e) d/A(R,p;£) (14a) 

/(R,p;e)= 2i5(p;e)6}A(R,p;e), (14b) 

where 

Of = [-2ze + v/(p)(V^-z^A(R))]-\ (15a) 

6} = h2ze~v/(p)(V^ + i— A(R))]"i. (15b) 



The strategy is to use a vortex lattice solution as an in- 
put, compute the anomalous Green's functions / and / in 
terms of g from Eq. (|14p . determine g from the normal- 
ization condition, and then enforce the self-consistency 
on A and the impurity self-energies. In principle, any 
complete set of basis functions is suitable for expanding 
both A(i?, p) and f{R,p). In practice, of course, we 
are looking for an expansion that can be truncated af- 
ter very few terms, enabling efficient computation of the 
functions. The Abrikosov lattice ansatz for A{R) is a su- 
perposition of the functions corresponding to the single 
vortex solution of the Ginzburg-Landau equations, and 
therefore it is natural to use these functions as our basis. 

For an s-wave superconductor with an axisymmetric 
Fermi surface (isotropic in the plane normal to the field), 
it is well known that the vortex lattice is given by a su- 
perposition of the single flux line solutions, the oscillator 
(Landau level, or LL) functions, ^o{x — xo), centered at 
different points in the plane normal to the applied field^ 

A{R)^Y.a^e^''«y^o (^^Z^^ . (16) 

ky 

Here the symmetry of the coefficients Cfc determines the 
structure of the lattice. This form emerges from the so- 
lution of the linearized Ginzburg-Landau (GL), and is 
also consistent with the solution of the linearized, with 
g = —in, quasiclassical equations. Moreover, this form is 
valid down to low fields as the admixture of the contribu- 
tions from higher Landau levels, $„ with n 7^ 0, to A(i?) 
remains negligible4^ Consequently, the set of oscillator 
functions, $„, provides a convenient basis for the expan- 
sion of anomalous functions /. It is common to rewrite 
the operator O via the bosonic creation and annihilation 
operators, a' and At the microscopic level, insert- 
ing this ansatz for A{R) into the quasiclassical equa- 
tions, Eqs. p4|) . and enforcing the self-consistency con- 
dition, yields the order parameter which only includes 
the ground state oscillator functions, justifying use of 
Eq. IHl)^. 

In unconventional superconductors the situation is 
more complex. While the solution of the GL equations 
arc still given by Eq. (fTB|) . this form is not a self- 
consistent solution of the linearized microscopic equa- 
tions: the momentum and the real space dependence of 
the order parameter are coupled via the action of the op- 
erator v/(p) • Vr,. Since the wave functions for Landau 
levels form a complete set, they can still be used as a 
basis for the expansion. The microscopic equations mix 
different Landau levels, and the self-consistent solution 
for the vortex state involves a linear combination of an 
infinite number of $„ at each site^. For the axisymmet- 
ric case the spatial structure of A(il) is still close to that 
for the s-wave case, and the weight of the higher Lan- 
dau levels in the self-consistent solution decreases rapidly 
with increasing ri^^^. Hence in practice the series in n 
is truncated either at n = (as for s-wave) or at the 
second non- vanishing ter m^^i^° . While this is often suffi- 
cient to describe the salient features of the thermal and 
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transport cocfRcients, care should be taken in determin- 
ing the anisotropies of these coefRcients under a rotated 
field: the anisotropy is often of the order of a percent, 
and the structure of the vortex lattice should therefore 
be determined to high accuracy as well. 

The situation is even more complex for unconventional 
superconductors with non-spherical Fermi surface, when 
the Fermi velocity is anisotropic in the plane normal to 
the applied field. Quasi-two dimensional systems with 
the field in the plane, such as shown in Fig.[l] give one ex- 
ample of such difficulties. Frequently in the microscopic 
theory the expansion is still carried out in the LL func- 
tions using the operators for the isotropic case. These 
functions are now strongly mixed, and hence (numeri- 
cally intensive) inclusion of many LL is required before 
the self-consistency is reached. Determining magnetiza- 
tion in the vortex state, for example, was carried out with 
6 LL functions^. 

This difficulty, however, is largely self-inflicted since, 
in contrast to the isotropic case, the LL functions in the 
form used in Ref.|4^are noi the solutions to the linearized 
GL equations. For an arbitrary Fermi surface the coeffi- 
cients of the Ginzburg-Landau expansion are anisotropic, 
and the vortex lattice solution is given by the n = Lan- 
dau Level in the rescaled, according to the anisotropy, 
coordinates4i We show in Appendix |^ that the proper 
rescaling is 



(17) 



lySf, y =y\/Sf 



where 5*/ is a measure of the anisotropy of the Fermi 
surface. For a FS with rotational symmetry around the 
axis zo, and for the field at an angle Oh to this axis. 



Sf — \ cos^ 6h H — 2^ sin^ 



Here 
;2 



and 



(18) 



''oil 



2(3^^(p)i^|(P2))fs, where v\\ is the projection of the 
Fermi velocity on the zq axis, and v^i with i ~ xo,yo 
is the projection on the axes in the plane normal to zq. 
For the field in the basal plane 9h = 7r/2, and therefore 

= voi\/vo±. 

The appropriate basis functions, which we use here- 
after, correspond to the oscillator states in the rescaled 
coordinates. If we chose the direction of the field as the 



z-axis, 



$„(a;,fcy) = 



A^/S 



(19) 



/ 



For an s-wave superconductor the n = ansatz for A(i2) 
satisfies microscopic equations, while for unconventional 
order parameters different LLs are once again mixed. 
However, with our choice of the basis functions this mix- 
ing is weak, enabling us to truncate the expansion at 
three components. Consequently, we use a generalized 



form of the vortex lattice A(il, p) = A{R)y{p), where 
A(i2) =^A„(R|??,) (20a) 

n 

(R|-> = E<^^7^'f"(-'^^)- (20b) 



The normalizing factor in Eq. (|20b[) is introduced so that 
the states ( R I n ) are orthonormal, i.e. 



dR 



Rl 



provided 



(21) 



(22) 



Consequently A„ in Eq. (j20ap has the meaning of the am- 
plitude of the appropriate component of the order param- 
eter in the LL expansion. 
The ladder operators, 

a = A (^_v,. + iiVy, + )) , (23a) 

at = A(^V.,+z(V,-+z|^)) , (23b) 

obey the usual bosonic commutation relations, [a, a^] = 
1, [a, a] [a^,a^] = 0, and connect the states |7i) via 
a\n) — ^/n\ n — 1 ) and a^\n) = ^/n + 1\ n-\-l). 

To solve Eg. p^ we rewrite the differential operators 
Of and Of via the ladder operators, Eq. (|^51) . and find 



Oi 



where 



-2ie-f v/(p) (V^-i^A(R))]- 



-2i£- 



V2A 



(u_(p)at - v+{i>)a) 



(24) 



(25) 



v±^v.,{t,)/^f±ivy{p)^f (26) 
For convenience we introduce the rescaled Fermi velocity 

Vf{-p)x=Vf{-p)^/^f , Vf(-p)y=Vf{-p)y^f, 

(27) 

and its projection on the I'y-plane (perpendicular to H), 

\vjm = ^vf{p)l+if{p)l, (28) 

as well as the "phase factors" , 

Vf{-p)x±ivf{v>)y 



■5±(p) 



(29) 



The off-diagonal parts of the matrix Green's function 
can be expressed in terms of the normal component 5, 
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and written as a series over the set ( R | m ) . The so- 
lution is based on exponentiating the operator O/ to 
exphcitly evaluate the result of its action on the order 
paramete r'^°i^^ , and is detailed in Appendix [Bl We find 



/(R,p;£) = ^/„,(p,e)(R| 



(30a) 



where w = \v^\/\/2A, and the prime over the k,l- 
sum denotes the restriction that the matrix element 
( n jaf^'a' | m ) = y^nlml/{n — k)l(m — l)\ is non-zero only 
for k < n, I < m and k — I = n — m. 

If we truncate the expansion of the order parameter 
in the vortex state at the lowest Landau level function, 
ri = 0; wc find from Eqs. ([33l) 



Uip,e) = ig ^(-«_(p)r-" |p|)A„(p;£) , 

n 

(30b) 

where A„(p;e) = A„(p) -I- H^imp.ni^)- The coefficients 

min[ra^n) 



2?m,n(£, |P 



(31) 



with ni(j) = 7 + (|m — n| — (m — n))/2, n2(j) = j + (I'ti- 
n| -t- {m — n))/2 in each term and 



^"'"'^^^"'.71/ (n-ni)!ni!n2 



(32) 

where VF("^(z) is the n-th derivative of the function 
W{z) = exp(— 2;^)erfc(— «z). These functions have the 
following symmetries: VF(")(z)* = (-l)"W^(")(-z*), 
= (-1)"-«I?„,„ and i?;^%"^(z)* = D:^\';:-{-z*). 
The diagonal part, g, is determined from the average 
// and the normalization condition. The details, once 
again, are relegated to Appendix [Bl with the result 



-in/Vl + P. 



(33a) 



x{n\a'^''a'-\m) 



n m k,l>0 
k+l 



l\k\ 



x/2 



1) / V2e 



, (33b) 



(34) 



1 - 



2A 



W^'(^)AoAo 



with previously obtained 



which agrees 

2 

expressions^ 

A cx oo, we use the asymptotic behavior 

at large values of the argument, W{z) « i/^/irz, 
W'{z) w — ^/(y^z^), to verify that this Green's function 
tends to the ECS limi t^^i^^ , and therefore all the 
conventional results for the density of states in nodal 
superconductors immediately follow. 

Eqs. ipO)) and ([55)1 give the solution of the quasiclassical 
equations in the BPT approximation for a given vortex 
lattice and impurity self-energies, i.e. provided the coef- 
ficients A„, A^ and T,, Aijnp,n, ^imp,m s-^c known. The 
self-consistency equations for these coefficients. 



A„ In ■ 



T 



TcQ J 47ri 



dpps nf{p)y{p) 



(35) 



2.^^^Vanh^, 
2T ' 



and the equations for the impurity retarded and ad- 
vanced self-energy, Eq. written explicitly through 
solution of Eq. (|10p for the i-matrix, 



r sin^ (5o 



"^p 1- ^((5)2 -(/)(/) +7r2) 



COtSo + {g)/TT {{f)/TT)i(J2 

«cr2((/)/7r) cot(5o - (5)/7r 



(36) 



complete the closed form solution. Here Tco is the critical 
temperature for the clean system, F = 0, which we used 
to eliminate the interaction strength, Vs, and the high 
energy cutoff. The elimination can also be done in favor 
of the impurity suppressed Tc, see e.g. Refl49l 



IV. HEAT CAPACITY 

A. Density of states and the specific lieat 

Once we self-consistently determined the Green's func- 
tion, we can calculate the quasiparticlc spectrum. We use 
the standard definition for the angle-resolved density of 
states at the Fermi surface, 



Nfip) 



1 



3mg^(p,e) , 



(37) 
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where Nf is the normal state DOS. 

The heat capacity is the derivative of the entropy, C 
TdS/dT, where 



S = -2^[(l-/(i?k))ln(l-/(i?k)) + /(i?k)ln/(i?k)] 

k 

-\-QQ 

= -2 [ deiV(e)[(l-/(e))ln(l-/(e)) + /(£)ln/(e)] 



/(e) = l/(e'^/"^ + 1) is the Fermi function, and N{e) = 
/ dp N{e, p) is the net DOS at energy e. In practice, nu- 
merical differentiation of the entropy is computationally 
either noisy or very time consuming due to high accuracy 
required in finding S, and therefore not very convenient. 
At low temperatures the order parameter and the density 
of states are weakly temperature dependent, and there- 
fore the specific heat can be obtained by differentiating 
only the Fermi functions. This leads to the well-known 
expression 



CiT,H) 



1 



-\-oo 



de 



e"^ N{T,H;e) 
cosh^(e/2r) 



(38) 



that lends itself more efficiently to numerical work. Note 
that the x^/ cosh^(a;/2) function has a single sharp peak 
at a; ~ 2.5, so the DOS at e 2.5 ~3T contributes the 
most to the C/T. The difference between the specific 
heat defined from the density of state and the exact re- 
sult is, of course, dramatic near the phase transition from 
the normal metal to a superconductor, where the peak 
in the specific heat is entirely due to entropy change not 
accounted for in Eq. ([55)) . At the same time, the regime 
where the anisotropy of C (T, H) is measured is far from 
Tc, and there we find that the results are very weakly de- 
pendent on the method of calculation. We therefore use 
the approximate expression above except where noted, 
and give a more detailed account of the difference be- 
tween the two approaches for the specific Fermi surface 
shape in Sec. IIV CI 



B. Cylindrical Fermi surface 

We are now prepared to consider the behavior of the 
specific heat in the vortex state of a superconductor. As 
mentioned above, our goal is to analyze the variations of 
the specific heat when the applied field is rotated with 
respect to the nodal directions. We consider first the 
simplest model of a cylindrical Fermi surface with vertical 
lines of nodes, and the field applied in the basal plane, 
at varying angle to the crystal axes. 

This is a simplified version of a model for layered com- 
pounds, such as CeCoIns, considered below in Sec. IIV Cl 



There we compute the specific heat for a quasi-cylindrical 
Fermi surface, open and modulated along the zg-axis. 
The main advantage of considering an uncorrugated 
cylinder first is that it provides a good basis for semi- 
analytical understanding of the main features of the ther- 
modynamic properties. Moreover, this model gives re- 
sults that are in semi-quantitative agreement with those 
for the more realistic model of Sec. IIV CI 

The disadvantage of the model is that it is not self- 
consistent. If the Fermi surface is cylindrical, there is 
no component of the quasiparticlc velocity along the zq 
direction (the axis of the cylinder). The field applied in 
the plane does not result in the Abrikosov vortex state, as 
the supercurrents cannot flow between the layers. Con- 
sequently, it is impossible to set up and solve the self- 
consistency equations for the order parameter as a func- 
tion of the applied field. Nonetheless we assume the exis- 
tence of the vortex lattice where the order parameter has 
a single n = Landau level component, with the ampli- 
tude MT,H) = A(r,0)v/1 - H/Hc2{T), analogous to 
Refiii. With this assumption, we solve self-consistently 
for the temperature-dependent A(T, 0), and for the im- 
purity self-energies. We consider the unitarity limit of 
impurity scattering (phase shift So = ""/S). In the next 
section we compare this model with a more realistic fully 
self-consistent approach, and show that the major fea- 
tures of the two are very similar. 

While in the cylindrical approximation the results de- 
pend solely on the ratio H/Hc2, for comparison with the 
results of the self-consistent calculation we recast them 
in similar form. We measure the field in the units of 
Bq = <I'o/27r^o where $o ~ hc/2\e\ is the flux quan- 
tum and ^0 = hvf/2'KTc is the temperature indepen- 
dent coherence length in the a6-plane. At zero tempera- 
ture the upper critical fleld along the c-axis is computed 
self-consistently, Hc2,c ~ 0.55iJo. We set the in-plane 
Hc2 = 1.1 -Bq to approximate the factor of 2 anisotropy 
found in CeCoIus, and choose the normal state scattering 
rate T /2'kTc = 0.007 (suppression of the critical temper- 
ature (Tco - Tc)/TcQ « 5%). We checked that the result- 
ing map of the anisotropy in the specific heat in the T-H 
plane does not strongly depend on this particular choice. 
Of course, large impurity scattering smears the angular 
variations. 

For a single Landau level component the solutions for 
the Green's function have a particularly simple form of 
Eq. (|34p . For a d^ja-ya superconductor the gap function 
is A(0) = Acos2</). If the magnetic field is applied at an 
angle 0o to the x-axis (inset in Fig. [2|) . the component of 
the Fermi velocity normal to the field is 

«^(0) = «/sin(0-(/.o). (39) 

Therefore the Green's fimction of Eq. ljM]) takes the form 
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9{e, ^) = (40) 



Let us focus first on tlie residual density of states, 
e ^ 0+ in the clean limit, F = 0, to compare with 
the semiclassical Doppler approximation. In this case 
W (0) = 2i/ and the density of states reduces to 



Heat capacity: C/y T 



7V(0) 



2tt 



d(f> 



1 



Q 2n |^ I 1 cos^ 2^ 



(41) 



where z = u//4\/2AA - ^H/Hc2- The DOS can be 
obtained analytically for the nodal and antinodal align- 
ments of the field. 

Node, 00 = 7i'/4. Then Eq.lgT]) reduces to 



K 



1 



(42) 



where K is the complete elliptic integral of the first kind. 
We use the convention of Ref . [10 for the argument of all 
elliptic functions. In the weak field limit, z <C 1, 



2? 4 
^„„,.(0)c.-ln-. 

TT Z 



(43) 



Antinode, (po = 0. The corresponding DOS is evalu- 
ated to be 



iVa„ti„odc(0) = 



(z2 + 1/4)1/4 TT 



K{r)-\F{a,r) 



(44a) 

where F{a^ r) is the incomplete elliptic integral of the 
first kind, and 



a 



1 - v/z^ + 1/4 




l + z2 



1 + + 1/4 ' \/2 Y ^ ' 
At low fields, ^ <C 1, the antinodal DOS 



= (0) 



2V2z, 4^2 
In 



4z2 
(44b) 



(45) 



Apart from the logarithmic correction (which is rapidly 
washed away by finite impurity scattering), the antin- 
odal DOS exceeds the nodal value by a factor \/2, in 
complete agreement with the Doppler approach^. As 
the field increases however, Eqs. (|42|) and (|44p predict a 
crossing point z* ~ 0.63 above which the residual nodal 
DOS becomes greater than A^ant nodo(O); this result was 
obtained numer ically in Refs. QlSl- 

With our choice of 
A(iJ) = A ■\/]~~h7^^c2 , the zero-temperature crossover 
point hes at H*/Hc2 ^ 0.6. 

Similar analytic expressions cannot be written for fi- 
nite energies and we evaluate the DOS and the specific 
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FIG. 2: (Color online) Left panel: The phase diagram for the 
anisotropy of the heat capacity for cylindrical Fermi surface. 
At low T and H (shaded area) the minimum in the heat capac- 
ity occurs when the field point in a nodal direction, 0o ~ 45. 
As T increases the minimum first evolves into a maximum, 
and then switches back to a minimum. The inversion of zero 
energy DOS is indicated by the dotted line. Right panel: evo- 
lution of the heat capacity anisotropy with temperature for 
H/Hc2 = 0.136 (circles in the left panel). Some curves are 
shifted vertically for clarity, their original values at 0o = 
are shown in boxes. 7jv is the Sommerfeld coefficient in the 
normal state. 



heat numerically, including the impurity effects. Re- 
sults for the anisotropy of the heat capacity are shown 
in Fig. [21 We present them in a form of a phase dia- 
gram (left panel) that shows the regions with the op- 
posite anisotropy. Shaded (white) areas correspond to 
the minimum (maximum) of C when H is along a node. 
Of course the node-antinode anisotropy disappears as the 
field — > 0. Since we are primarily interested in compar- 
ison of our results with the experimental data, we focus 
on the regime of moderate fields, and show the evolution 
of specific heat for different directions of the field, 
with the temperature in the right panel of Fig. [2l Notice 
that at 00 = 45°, when the field is along a nodal direc- 
tion, the minimum in C(0o) evolves into the maximum 
as T increases. 

Inversion of the anisotropy in the T-H phase diagram 
is at odds with the semiclassical result that always pre- 
dicts a minimum in the specific heat for the field parallel 
to the nodal direction. In the shaded area adjacent to the 
Hc2(T) line in Fig. [2l with minima for H|| node, the spe- 
cific heat is already sensitive to the density of states near 
the BCS singularity in the DOS at e '--^ Ao, and there- 
fore direct comparison with the semiclassical analysis is 
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N(e; H) / N 




FIG. 3: (Color online) Left: the low-energy part of the DOS 
for cylindrical FS. The nodal and antinodal DOS cross at 
finite energies (arrows). Right: the angle-resolved DOS (red 
shaded) for the two field orientations in the regions indicated 
by the dotted boxes in the left panel. The angle integrated 
DOS is given by the area of the shaded regions. See text for 
details. 



not possible. Moreover, we show in the following section 
that the self-consistent models require nodal-antinodal 
anisotropy of the upper critical field, and the results for 
this part of the phase diagram are modified. 

On the other hand, the anisotropy inversion between 
the low-T, low-H region, and the intermediate temper- 
atures and fields, occurs still in the regime where the 
semiclassical logic may have been expected to work. The 
dotted line in the left panel of Fig. [5] separates the two 
regions of the residual zero-energy DOS: below that line 
iV„ode(0) < iV,„ti„odc(0), while above the line iV„odc(0) > 
-^antinodo(O). Thc iuvcrsion of the anisotropy in the spe- 
cific heat is clearly not just a consequence of the behavior 
of the zero-energy DOS found above. Recalling that C/T 
is predominantly sensitive to the density of states at en- 
ergies of the order of a few times T (see Eq. ([55|) ). we 
conclude that the origin of the anisotropy inversion is in 
the behavior of the finite energy DOS. We plot the low- 
energy N{e, H) at several values of the magnetic field in 
the left panel of Fig. [3] At low fields, the DOS anisotropy 
at small e agrees with the semiclassical prediction, but 
the density of states for the field along a node (dashed 
lines) and along an antinodc (solid lines) become equal 
at a finite energy indicated by arrows. Above this energy 
the DOS anisotropy is reversed, and is manifested in the 
reversal of the specific heat anisotropy as T increases. 
The crossing point moves to lower energies with increas- 
ing field, and is driven to zero when the residual, e = 0, 
DOS for the two directions become equal. In our numer- 
ical work with finite impurity scattering rate this occurs 
at H* 0.5Hc2, and we checked that H* 0.6i7c2 as 
the system becomes more pure, in agreement with the 
analytical results above. 

As suggested by us in the short Letter communicat- 



ing our main results, the inversion stems from the in- 
terplay between the energy shift and scattering due to 
magnetic field^. Magnetic field not only creates new 
quasiparticle states on the Fermi surface, but also scat- 
ters the quasiparticles and, consequently, re-distributes 
their spectral density. This scattering is present in 
the microscopic method, but not in the Doppler shift 
treatment. To understand this effect and to make con- 
nection with the semiclassical approach we analyze the 
angle-resolved DOS obtained from the Green's function, 
Eq. (p4|) . It is instructive to re- write the Green's function 
in the BCS-like form which makes explicit the distinc- 
tion between the energy shift and scattering rate. We 
define the "magnetic self-energy" = — iE" from 
{e- E)-2 = i^{2A/\df\)^W'{2iA/\ijf\) so that the 
Green's function reads 



|Aop3^2(p) 



(£-I]'(e,H,p)+zE"(£,H,p))2 



-1/2 



(46) 

The density of states for a given direction at the Fermi 
surface can be found from the comparison of £ — E(p) 
with Ao3^(p) = Amax cos 20. Since W'{x) is a complex- 
valued function, both E' and S" are generally non-zero: 
the former shifts the quasiparticle energy, while the lat- 
ter accounts for the direction-dependent scattering. For 
now we neglect the impurity broadening: for quasiparti- 
cles moving not too close to the field direction the field- 
induced scattering is normally greater than the scatter- 
ing by impurities. Non-zero E" is the key signature of 
our microscopic solution. Both real and imaginary com- 
ponents of the self-energy depend on the quasiparticle 
energy e, the strength and direction of the field H and 
on the momentum of the quasiparticle with respect to 
both nodal direction and the field. Using the expansion 
around W'{0) = 2i/y/n at small values of the argument, 
and taking W'{z ^ 1) « —i/y/nz'^ for large arguments, 
we find two limiting cases, 

HI 



e-E 




if£»^ - (48) 



H is the characteristic mag- 



Note that 

netic energy scale for quasiparticles at position p on the 
Fermi surface. In the first limit, valid at low energies 
(or moderately strong fields) for quasiparticle momenta 
away from the field direction, the imaginary part of the 
self-energy is dominant. In the opposite limit the effect 
of the field is small. Between these two limits, i.e. at 
finite energies, moderate fields and arbitrary p , the real 
(energy shift) and the imaginary (scattering) parts of the 
self-energy can be comparable. 

We can now analyse the angle-dependent contribution 
to the density of states from different regions at the Fermi 
surface at a given field, which is shown in the right panel 
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of Fig. [3l Consider first very low energy e ^ 0, panels 
Fig. [21b) and c), so that we are in the regime described by 
Eq. (j47ll . At low fields, panel c), the characteristic energy, 
|u^(p)|/2A is smaller than the maximal gap, Amax, and 
therefore most of the field-induced quasiparticle states 
appear near the nodes for which [^^(p)! is moderately 
large. Consequently, as in the semiclassical result, field 
applied along a nodal direction does not create quasipar- 
ticles near that node, while the field applied along the 
gap maximum generates new states at all nodes. The 
small contribution seen in the right frame of panel c) at 
the nodes aligned with the field is due to impurity scat- 
tering. Thus, while the scattering on the vortices, i.e. the 
imaginary part of Eq. ^47\i , does produce a non- vanishing 
contribution to the DOS over most of the Fermi surface, 
at very low energy and low field the spectral weight of 
the field-induced states is mainly concentrated near the 
nodal points. 

This changes as the field is increased, see panel b). 
At high field the Doppler shift and pairbreaking due to 
scattering are strong, and sufficient to contribute to the 
single particle DOS over almost the entire Fermi surface 
where |'u^(p)|/2A ^ Amax{H) (as a reminder, in our no- 
tations the maximal gap, Amax = Aq^/2, since we chose 
yicp) = cos 20 to be normahzed). The obvious excep- 
tions are the momenta close to the direction of the field, 
when |w^(p)| ^ Vf. For the field aligned with the node 

this restriction is not severe: near the node ~ VfS(l) 
and A(p) ~ 2AmaxS(l), where Scj) is the deviation from the 
nodal (and field) direction. Hence if Vf/{2A) ~ A^ax, 
almost the entire Fermi surface contributes to the DOS 
when the field is aligned with a node. In contrast, for the 
field along the gap maximum, in a range of angles close 
to the field direction the gap is large and the magnetic 
self-energy is small, and hence no spectral weight is gen- 
erated. As a result, the density of states is higher for the 
nodal orientation. This is the origin of zero-energy DOS 
inversion as found numerically in Ref. H and as derived 
above. 

We finally consider the DOS at finite energy e ^ 
Amax, panel (d). In the absence of the field the most 
significant contribution to N{e) comes from the BCS 
peaks at e = Ao|3^(p)|, located at momenta p^ at an- 
gles (f)n ± Scj);,, where 0„ = 7r/4 + Tm/2 are the nodal 
angles, and (50e w e/{2Afnax) = e/(2\/2Ao). Scatter- 
ing on impurities or vortices broadens these peaks, and 
re-distributes the spectral density to different energies 
(as in all unconventional superconductors, scattering re- 
duces the weight of the singularity and piles up spectral 
weight at low energies). However, the vortex scattering 
is anisotropic as it depends on v^, see Eq. (j47| . the com- 
ponent of the velocity normal to the field. Therefore if 
a field is applied along a nodal direction, at that node 
vj- ~ Vfd(j)e <^ Vf, and the peaks in the angle resolved 
DOS remain largely intact (d, right). On the other hand, 
if the field is applied along a gap maximum, BCS peaks 
near all four nodes are broadened by scattering, and their 



contribution to the net DOS is reduced (d,left). So, even 
when the field is moderately low but the quasiparticle 
energy exceeds some value e*, which can only be de- 
termined numerically, the gain from sharp (unbroadened 
by scattering) coherence peaks exceeds the field-induced 
contribution from the near-nodal regions. Then the DOS 
is higher for the field along a node rather than the gap 
maximum. Recalling that the specific heat at tempera- 
ture T is largely controlled by the density of states at the 
energy of about 2.5T, we expect that the anisotropy of 
the specific heat is also inverted at T/Tc ^ e* /2.5Tc. It 
is this change in the finite-energy density of states, rather 
than the zero energy DOS, that determines the inversion 
line in the phase diagram, see Fig. [2] 

C. Quasi-two-dimensional Fermi surface 

We mentioned above that a major motivation of our 
work is to address the apparent discrepancy between the 
thermal conductivity and specific heat measurements in 
CeCoIns. While this material does possess a quasi- two 
dimensional sheet of the Fermi surface, the normal state 
resistivity anisotropy is very moderate, indicating a sig- 
nificant c-axis electronic dispersion. Consequently, while 
the results of the previous section are very suggestive of 
the anisotropy reversal, we need to verify that similar 
physics persists in a more realistic open quasi-cylindrical 
Fermi surface, described by 

P} "^^pI+pI^ {r'^ p))cos{2spJr'^Pf) . 

We parameterize this FS by the azimuthal angle in the 
a6-plane, 0, and momentum along the c-axis, Pz, so that 
the Fermi velocity at a point ((j),pz) is 

(Pf/m + cos{2. spz/r'^Pf) costj) \ 
Pf/m ^1 + cos(2spz/r^Pf) sin0 
Pfs/m sm{2spz/r'^Pf) ) 

With this parametrization, the anisotropy factor of the 
normal state DOS is nf{p) = 1. Parameter r deter- 
mines the corrugation amplitude along the z-axis, and 
we find that the results do not depend on its value; be- 
low we set r = 0.5. The second parameter, s, is physi- 
cally important since it fixes the anisotropies of the nor- 
mal state transport and the critical field: the character- 
istic velocities in the a5-plane and along the c-axis are 
''^O-L = Pf/iTT- = Vf and fo|| = Pfs/m. The normal state 
conductivity anisotropy is (Jzz/<^xx = s^- 

The main advantage of allowing for the c-axis disper- 
sion is the ability to solve the quasiclassical equations 
in the BPT approximation self-consistently, with respect 
to both the order parameter as a function of T, H, and 
the impurity self-energy. We take moderate values of the 
anisotropy, s = 0.25 and s = 0.5, for which the vor- 
tex structure is still three-dimensional. The latter value 
yields Hc2 anisotropy close to that of CeCoIs. The calcu- 
lations below are done with three Landau level channels 
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FIG. 4: (Color online) Multiple Landau level contributions 
to the order parameter. Left panel; LL components as a 
function of the field in the antinodal direction for T/Tc — 
0.2 (left). Right panel: upper critical field along a node for 
difi'erent number of Landau channels in A (clean limit); Hc2 
has converged for A'' > 3. We also show for comparison the 
Hc2 as calculated with one channel Ao without coordinate 
rescaling (|17fl . 



for the order parameter, Aq, A2, A4. With the rescaling 
of Appendix |X] this is sufficient for convergence of the 
upper critical field. The values of the higher components 
A2, A4 are less than 5% of Aq, see Fig. |4l and addition 
of further components does not change the results. 

For this Fermi surface we solve the linearized self- 
consistency equation and compute Hc2 in the basal plane. 
The anisotropy between nodal and antinodal upper criti- 
cal fields appears naturally as a result of the d-wave sym- 
metry, H^^'^" ^ H^ntinode^ r^^iQ valuc of i?c2 is essentially 
determined by balancing the kinetic energy of the super- 
currents vs. the condensation energy, and the former is 
different for different orientations of the field. 

Let us now look at the difference between the self- 
consistent and non-self-consistent order parameter cal- 
culations. For this we again present a phase diagram, 
Fig. [51 analogous to Fig[2] for the cylindrical FS. Left 
panel shows the results for the Fermi surface with r = 
s = 0.5, and the impurity strength V/2t:Tc = 0.007 
(Tc/Tco ~ 0.95, ^ 70). The values of the crit- 

ical fields at T = are i/™*™"'^'^ « 1.45So, H^i'^" « 
1.27i?o and H^2 ~ 0.57i3o- This gives the in-plane 

anisotropy [H^^^^'^^de _ jjnodeyjjantinode ^ ^^^^ 

the ratio between the c-axis and antinodal directions, 
■f^c2/-f^c2"*"°'^'' = 0.4. To demonstrate the influence 
of the FS c-axis curvature on this phase diagram, we 
present in Fig. [^b) a similar diagram a Fermi sur- 
face with parameters, r = 0.5 and s = 0.25. These 
parameters correspond to the reduction by factor of 2 
of the velocity along the c-axis and the critical fields 
jjanuuode _ 2.85Bo, H^.^'^'' « 2.55Bo and H^^ « 0.57Bo- 
The Hc2 anisotropies are: 10% in the basal plane between 
nodal and antinodal directions, and H^2/ ^^2*"''°'^'^ ~ 0.2. 
Fig. [5] shows that a factor of two difference in the c-axis 




T/T T/T 

c c 

FIG. 5: (Color online) The phase diagram of the heat capacity 
anisotropy C{(f>o) for the corrugated FS with s = r = 0.5 (left) 
and s = r/2 = 0.25 (right). 

velocity affects only the absolute values of ^^°"**)"°'*'^^ 
but otherwise the two diagrams for the anisotropy in the 
a6-plane look almost identical. 

The shaded "semiclassical" region at low temperatures 
and fields in Fig.O where minima of C are for H|| node, 
expanded compared with that for cylindrical FS (Fig. [2]). 
We note that if we truncate the order parameter ex- 
pansion at the lowest Landau level, without full conver- 
gence of Hc2, the "nodal minimum" region occupies sim- 
ilar range for both corrugated and purely cylindrical FS. 
Therefore this expanded range is the result of the self- 
consistency and inclusion of higher harmonics. On the 
other hand, the shaded 'minimum-at-a-node' region near 
Hc2 shrunk to low H and high T, where the anisotropy is 
almost washed out, and is experimentally undetectable. 

Specific heat as a function of the field direction is 
shown in Fig. [6l The curves are computed at the (T, H) 
points indicated in the phase diagram of Fig. [S] by cir- 
cles and squares. The left (right) panel refers to lower 
(higher) field. At higher fields the gap nodes always cor- 
respond to maxima of C. At low fields, however, nodes 
correspond to either minima or maxima of C depend- 
ing on the temperature. Fig. [B^left). The lowest dashed 
curve in the left panel of Fig. [S] appears to contradict the 
semiclassical results; we show below that it corresponds 
to the breakdown of the BPT approximation. We con- 
clude that the optimal range of field and temperature for 
experimental detection of the nodes based on the heat 
capacity anisotropy is at intermediate values of H/Hc2 
and T/Tc, where the anisotropy of C is large and the 
ambiguity in interpretation is small. 

The discrepancy between t = 0.5 profile on the left 
that shows a weak minimum at the nodes and the posi- 
tion of the point in the 'maximum' region of the phase 
diagram in Fig. [5] is due to the fact that we computed 
and differentiated entropy to determine the phase dia- 
gram, but employed the approximate formula Eq. (|38p to 
calculate the heat capacity anisotropy profile (neglecting 
the derivative of the DOS with temperature). Compari- 
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FIG. 6: (Color online) Heat capacity profiles for T-scans at 
H= and H=. The gap node is at (jio = 45°. Left: the 
evolution of anisotropy with temperature along a line of circles 
(low field) in Fig[5ja). Right: same for the high field (squares 
in FigO^a)). The curves that were shifted down to appear on 
the same scale are shown with their original values at 00 = 
in boxes. 




T/T 

C 



FIG. 7: (Color online) The heat capacity at H/Hc2 = 0.27. 
Comparison of the approximate formula Eq. (|38p with the 
rigorous calculation of C/T = dS/dT from numerical differ- 
entiation of the entropy. The crossing of C's for nodal and 
antinodal H directions at low temperature T/Tc « 0.15, in- 
dicated by arrows in the inset, is not significantly affected by 
the approximation. 



son of the exact and approximate formulas for the heat 
capacity is shown in Fig. [71 The lower inversion between 
the minimum and the maximum of C for the field along 
the nodes is only slightly shifted to higher T due to use of 
the approximate formula (inset). The point of the high-T 
inversion is more sensitive to it, but, as discussed above, 
is not in the regime of experimental interest. 

At the lowest fields and temperatures in Fig. [5] (below 
Q.lHc2 and O.OTTc) there appears a very small anoma- 
lous region where the heat capacity anisotropy is inverted 
compared to the semiclassical result. Our analysis shows 
that this is an artefact caused by the breakdown of the 
BPT approximation. Manifestations of this failure are 




I , \ , \ , \ , I I , I , \ 1 \ , I L_i \ , \ , \ , I 
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E / 27tT^ E / 27lT^ E / 27tT^ 



FIG. 8: (Color online) The DOS for quasi-cylindrical FS 
for fields H/Bo = 0.1,0.4,0.8 {H^2 = 1.45B(,(1.27Bo) for 
H|| antinode (node)), (a) Overall shape of the DOS for three 
fields, and (b) enlarged low-energy region. The low field DOS 
for nodal (dashed) and antinodal (solid) directions of H cross 
several times. We ascribe it to high derivatives of H^-function. 
The oscillating behavior decreases as system becomes cleaner 
(c), and disappears for single channel (d). 



enhanced (compared to cylindrical FS) by the fully self- 
consistent calculation of the multiple Landau channel or- 
der parameter. 

A necessary condition for the validity of the BPT ap- 
proximation is that the electron mean free path is much 
greater than the intervortex distance, t{H) = VfT{H) ^ 
A(i/). Only in this case we are allowed to carry out 
the vortex lattice spatial averaging before averaging over 
the impurity configurations to compute the self-energy. 
Consequently, for finite impurity concentration the ap- 
proximation is bound to fail at low fields. In figure [5] we 
present the DOS at low field and temperature for differ- 
ent number of A-channels and the purity of the material. 
Notice that for the dirtier material with more than one 
channel of the order parameter the DOS oscillates at low 
energies when the field H is || antinode, panel (b). These 
oscillations lead to the additional unphysical inversion of 
the heat capacity anisotropy at very low T,H, that is 
seen in the bottom left corner of the phase diagram in 
Fig. [51 and is shown by the dashed line in the left panel 
of Fig. [HI The same oscillations are also present in the 
sclf-consistcntly calculated impurity self-energies, which 
we do not show here. We find that they decrease in a 
cleaner material. Fig. [5Jc) . 

The interval of the fields where the oscillations are 
observed coincides with the region where the BPT ap- 
proximation is no longer trustable. We consider the 
BPT breakdown at low temperature, where the den- 
sity of states is, Eq. ((301), N{0,H)/No = -3m{g/n) ~ 
■\/ H / Hc2- The impurity self-consistency in the unitary 

limit gives t{H) = ^ N{0,H)/No ~ ^ ^/^- Here 

7 = r sin^ So is the nor mal sta te scattering rate. Recall- 
ing that A{H)/^o ~ \flld2fii and requiring VfT{H) ^ 
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A{H), we obtain a condition for the applicability of BPT: 
H/Hc2 ^ 7/27rTc. Thus, for our impurity bandwidth 
^ /2'kTc ^ 0.01 the BPT approximation is only applica- 
ble for fields H/Hc2 S> 0.01 and the oscillations seen in 
the DOS likely are a signature of this breakdown. We 
checked that increasing disorder expands the anomalous 
region and is consistent with this interpretation. For the 
single-component (lowest Landau level) A, the numeri- 
cally computed DOS does not show significant anoma- 
lous behavior, Fig. [SJd), at the same impurity level as 
in Fig. [5Jb). We argue that although the breakdown 
of the approximation is still there, its manifestation is 
less pronounced compared with the multiple-channel or- 
der parameter. Use of higher Landau channels for the 
expansion of A leads to the appearance of the higher 
derivatives of the l/l^(2;)-function, in the Green's 

function, see fEa l33[) . These grow very fast with n at 
z = Q (approximately as n\\) and are strongly oscillating 
as the argument is increased from zero. This is likely the 
underlying reason for the oscillations in the DOS, and 
therefore the ultra-low T-H inversion is an artefact of 
using the approximation beyond its region of validity. In 
contrast, other inversion lines in the phase diagram cor- 
respond to physical inversion of the measured properties. 



V. DISCUSSION AND CONCLUSIONS 

In this work we laid the foundations for an approach 
that provides a highly flexible basis to the calculation, on 
equal footing, of the transport and thermodynamic prop- 
erties of unconventional superconductors under magnetic 
field. The theoretical method is based on the quasiclas- 
sical theory of superconductivity and the Brandt-Pesch- 
Tewordt approximation for treatment of the vortex state 
in superconductors. This approximation allows for ac- 
curate and straightforward (analytic closed form expres- 
sions for the Green's functions) way to describe effects 
of the magnetic field in almost the entire T-H phase 
diagram for clean superconductors, with the exception 
of ultralow fields and temperatures. Combined with the 
non-equilibrium Keldysh formulation of the quasiclassical 
theory it paves a path for a very effective computational 
scheme that self-consistently takes into account multi- 
ple Landau levels of the expansion of the order parame- 
ter and impurities, and allows calculations for arbitrary 
temperature and magnitude of the field. The compan- 
ion paper II extends the method to the calculation of 
transport properties and focuses on the electronic ther- 
mal conductivity. 

Here we computed the density of states and the specific 
heat in the T-H plane for a d-wave superconductor with a 
quasi-two dimensional Fermi surface (cylinder modulated 
along the symmetry axis), and the magnetic field rotated 
in the basal plane. This choice of the Fermi surface and 
the field orientation was motivated by experiments on the 
heavy ferniion CeCoIng^. We provided the first com- 
plete description of the evolution of the anisotropy of the 




FIG. 9: (Color online) The anisotropy of the heat capacity 
in rotated magnetic field. The T-H phase diagram for a d- 
wave superconductor with a corrugated cylinder Fermi sur- 
face with purely orbital depairing. In the large part of the 
phase diagram the maxima of C{4)o) function correspond the 
nodal directions (indicated by arrows) on the Fermi surface. 
This result is opposite to that in the Doppler region, which is 
confined to H/Hc2 < 0.5 and T/Tc < 0.2. 



heat capacity due to nodes of the superconducting gap 
across the T-H phase diagram, see Fig. [9l 

Our main conclusion is that the anisotropic scattering 
of quasiparticles due to vortices plays a crucial role in the 
variation of the density of states and the specific heat as 
a function of the field direction. This effect is absent in 
the semiclassical (Doppler shift) approach, and becomes 
important already at moderately low fields, and at finite 
temperatures. As our phase diagram of Fig. [9] shows, as 
a result of this scattering, the anisotropy in the specific 
heat changes sign as a function of T and H. At low fields 
and temperatures the minima in the heat capacity occur 
when the field is oriented along the nodal directions, in 
agreement with the semiclassical (Doppler shift) calcula- 
tion. At higher T and H (already at T/T^ > 0.2 at low 
fields) the situation is reversed, and the maxima rather 
than minima of the specific heat are found when the field 
is along a nodal direction. Moreover, we showed that 
the inversion is related to the behavior of the density of 
states at finite energy, and not simply the residual DOS 
at the Fermi surface. 

While we expect that the loci of the inversion lines in 
the T-H plane weakly depend on the shape of the Fermi 
surface^, it is the existence of this inversion and its con- 
nection to the scattering and the finite energy DOS that 
emerged from our theoretical description and was not 
captured by previous approaches to the problem. Our 
calculations serve as a basis for the analysis of the exper- 
imental data, and we note that the interpretation of ex- 
periments based on the low-field expectations of minima 
for the field along the nodes can lead to diametrically op- 
posite conclusions regarding the gap symmetry, depend- 
ing on the values of the field and temperature where the 
anisotropy is been measured. The results suggest that 
the amplitude of the anisotropy is the greatest at inter- 
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mediate temperatures and fields, and that it is desirable 
not only to measure the C anisotropy at a few tempera- 
tures and fields, but also determine its evolution over the 
phase diagram. 

As an example, we consider the data for CcCoIns from 
Ref. [13, and plot in Fig. [9] the points where the published 
data were taken. The measured C{4>o) shows minima for 
the field along the [100] and [010] directions {(j>o = 7rn/2 
with n = 0,1,2,3), at all three locations, with vanish- 
ingly small anisotropy at point 3. Points 2 and 3 are 
clearly in the region where maxima of C{(f)o) determine 
the nodes, and thus firmly point towards d^^-y^ sym- 
metry. Point 3 is also close to the inversion line, which 
explains small amplitude of the oscillations. Point 1, in 
contrast, is in the "semiclassical" region, and therefore 
the minima of C((/)o) for the field along the crystal axes 
may be more suggestive of a dxy symmetry. Wc note, 
however, that the exact location of the inversion line is 
sensitive to the exact shape of the Fermi surface, and 
changes between the calculations restricted to the lowest 
Landau level for cylindrical Fermi surface. Fig. [5] and the 
multicomponent quasi-2D case. Fig. O We argue there- 
fore that points 2 and 3 are more reliable indicators of 
the gap symmetry, and the results are more suggestive of 
the dx2-y2 gap. While such a conclusion purely from the 
specific heat data is not foolproof, we show in II that the 
dx2-y2 symmetry is also supported by the analysis of the 
heat transport anisotropy of Rcf . [lol . 

The microscopic approach, by its very nature, couples 
the gap symmetry with the shape of the Fermi surface. 
For that reason direct comparison of our results with 
other experimental data, for example in the borocarbides 
YNisBzGi^ and LuNiaBaGi^, is not possible. These sys- 
tems are essentially three dimensional, and the Fermi sur- 
face has no quasi-cylindrical sheets. Moreover, it is very 
likely that there is substantial gap modulation along the 
z-axis, and comparison should be made with both point 
and line node model a^^'^^ . While the argument for the 
change in the anisotropy due to scattering on the vortices 
is quite general, the position of the anisotropy inversion 
lines in the phase diagram (if any) is undoubtedly differ- 
ent from that found for the quasi-2D system, and such 
differences are known to occur in the zero-energy DOS—. 
Therefore we will consider the nodal structures of these 
systems separately in near future. 

To reiterate, the approach described in this work 
presents a powerful tool to study the gap symmetry in 
the unconventional superconductors taking into account 
their realistic Fermi surfaces. Our results serve as a basis 
for interpretation of experimental data, pointing towards 
a resolution of the discrepancy between the results of the 
specific heat and thermal conductivity measurements in 
CeCoIns, which is also addressed in II. The method de- 
veloped here can be easily generalized to include other 
Fermi surfaces, paramagnetic effects, and other aspects 
of real materials, the discussion of which we defer to fu- 
ture publications. 
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APPENDIX A: CHOICE OF OPERATORS FOR 
ANISOTROPIC FERMI SURFACE 

The raising and lowering operators for the eigenfunc- 
tion expansion of the order parameter, a\ a and the cor- 
responding ladder states can be introduced in several dif- 
ferent ways. We want to define them in the manner that 
facilitates the efficient computations. This issue becomes 
important for anisotropic Fermi surfaces and arbitrary 
direction of the field. Anisotropy of the FS is directly 
translated into the shape of a single vortex and we can 
choose the orthogonal states such that they approximate 
this shape well already at the lowest order truncation of 
the expansion of A(i2). 

We consider an axisymmetric FS in cylindrical co- 
ordinates (r, (/), zq). The energy has the form 2ms — 
Pr + f{Pzo), with an arbitrary function fipzo) of the Pza 
momentum. The vortex state near the critical tempera- 
ture Tc is determined from the linearized GL equations, 

- i^i,- V.V, A(R) + (|r - l) A(R) = (Al) 

V = V - i— A(R) (A3) 
c 

In these equations the coordinates {x, y, z) are chosen 
so that the field is along the z-axis, B = Bi, and we 
take A = (0, i3a;,0). The form of the Kij tensor de- 
pends on the shape of the Fermi surface, the pairing 
state and orientation of the magnetic field. If the rota- 
tional symmetry axis is zq j the velocity of quasiparticles is 

V/(</',Pzo) = {vT{Pzo)cOS(l>,Vr{j)zo)^va.<P,Vzo(Pzo))- If B is 

along one of the FS symmetry axes, xq, yo or zq, Kij is di- 
agonal for d-wave pairing with y = \/2cos2(0 — 0o). We 
have, Kxgx„ = Kygy„ = Kov^^, Kzgza = -^oWoH' where 
Ko = 7C(3)/8(2^re)2, and 

vl^ - 2{y\i>)vl^^y^^){pz,)U, (A4) 
^^011 - 2{y\f>)vl{pz,)U. (A5) 

We apply the magnetic field at a tilt angle Oh from 
zq direction towards Xq axis. A coordinate system as- 
sociated with B is chosen as follows, z is along B, 
y = yo and x lies in (xq, zo)-plane and perpendicular 
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FIG. 10: (Color online) Quasi-cylindrical Fermi surface con- 
sidered in this paper. Direction of H defines ixyz) coordinates 
with z||H. To go from (xoyozo) coordiantes, associated with 
the Fermi surface, to {xyz) coordinates, associated with the 
field, we perform first rotation by (/)o around zq, and then 
rotation by around y. 



to z. Projections of the Fermi velocity at different points 
of the FS on these new coordinate axes, Vf^xiPzoi^') = 

zot4') — 1^f,zoiPzoT 4>) '^Os9h + 

Vf^xoiPzoj4>)siT3.0H- In (x, ?/, z)-coordinates the tensor 
Kij is not diagonal anymore, Kyy 



Kxoxo cos^ 9h 
K^„^„ cos'^ 9h, K. 



, sm tlH, 



6h. iC. 



Kxoxo sin^ 9h 



— (Kxgxo ~ KzQZQ)sai9H cosOh, and 
for the choice of the operator, V — z2e/c A = (Vj^, Vy — 
i2e/cBx, Vz), the GL equation is 



KxxVi A 



yy 



-2KxzVxV, A 



.2eS , ^ 
I X A 



— -1 1 A = 0. 

Tr 



(A6) 



It is easy to check by setting A ~ A{x,y) exp{ikzz) 
that the highest critical field still corresponds to fcz = 0, 
and we put V^A = below. We rescale the coordinates 
x' = x/^y^ and y' = j/y^^, and choose the scaling 
factor Sf such that Kxx/Sf = KyySf. Thus, 

,,2 



ri2 ^xx 

^^yy 



V, 



cos Oh 



V, 



^sin^ 9h. 

0_L 



(A7) 



After introducing creation and annihilation operators 
(A2 = c/2\e\B and e < 0), 



equation 



6]) becomes 



^(^y'+^j2 



(A8) 
(A9) 



A{x,y) 



A2 S 



f 



2 



1-^] A{x,y). 



(AlO) 



Then for any axisymmetric FS we obtain the well-known 
result of the anisotropic mass model for the upper critical 
field, determined by the ratio of the Fermi velocities for 
the two directions. 



Bc2{9h,T) 



const ■ (1 — T/Tc) 



cos2 9h + ^ sin^ 0H 



(All) 



We also find a set of eigenfunctions. 



(A12) 

In terms of the operators a, the gradient term in the 
Eilcnberger equation has the form 

2e 



v/(p) (V,-^-A(R)) = 
-^[-{)+(p)a + {)_(p) a}\ 



(A13) 



Here we rescaled the Fermi velocity in the zy-plane 



if{p)y = Vf{p)yy^, 



with 



and 



I^/(P)I = ^Jifip)l+Vf{p)l., 
ifip)x±ivfip)y 



v±ip) 



\vf\ 



(A14) 
(A15) 



(A16) 



(A17) 



APPENDIX B: CLOSED FORM SOLUTION FOR 
THE GREEN'S FUNCTION 



To solve the scmiclassical equations wc use Eq. (|A13[) 
to cast the operator O f from Eq. (|15p in an integral form, 



^ V2A 



where we introduced the magnetic field energy 



w 



V2A 



(B2) 
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and used the operator identity exp(j4 + B) = 
exp(A) exp(i3) exp(— C/2), if [A,B] = C is a c-number. 
In all integrals we also keep in mind that 2fm£ > for 
retarded functions, so that the convergence is ensured. 

In this formulation it is convenient to work with bra- 
and ket-functions for different vortex states, which cor- 
respond in R- representation to states (jA12[) . We present 
decomposition of A as 



(B3) 



with operator Of acting to the right, 
/ = {2tg)dfA = 



(B4) 



{2ig) J di^e2*^"*i-'"'*i/2e-"'*i"-''^e"'*i"+"A, 



and operator 0| acting to the left, 



(B5) 



/=(2zg)A6}. 
We rewrite the operator Of as 

Of ^ [2ie* -w{v-a^ -i+a)]-^ ^6^ ^ (B6) 

/ (ft g — 2iE*t2— 10^*2/2^10*2 S-O* g — '!«t2 -iJ+a 





so that the spatial average of the off-diagonal functions 
is 

00 00 

77 = {2igf J dh J dt2 e2.e(ti+t2)-» = (ti+t2)V2 X 


X ^ ^g— !i'(tl+t2)'0-o''"gio(tl+t2)'0+a~j ^ ^ (-gy^ 

where we make sure that bra-vectors stay on the left 
of ket-vectors. Here we again used operator-in-cxponent 
rule to commute exponents. After an appropriate vari- 
able susbstitution, 

00 

Yf = {2igf J dUe2*"*-™'*'/2A (^6-""*"-"^'"*"+°) A. 

(B8) 

This form is very convenient if we intend to keep several 
Landau channels in the expansion of A. If the highest 
Landau level used is N, the series expansion for 6™* "+" 
contains only N + 1 terms; and to calculate the spatial 
average / / we need to compute only a finite number, 
{2N + 1), of W^(")-functions, since 



dt t (wt)" e 



^1) / V2e 



(B9) 



Solution for / is written as, 

/(R,p;e) = E/™(p,£)(R|m), (BIO) 



with the amplitudes 
fmiP, e) - ig E(-«_)"-"(p) P„,„(£, |p|)A„(p; e) . 

n 

(Bll) 



Here 



25m,n(e, IpI) = 



D 



ni,n2/ .,\ _ 



(^) 



71 



ni+n2 



ni +712)! 
{n — 7ii)!7ii!n2! 



(B12) 



(B13) 



The sum starts from 712 = max(0, m — n) and, in each 
term, ni = n — m -\- 712- This sum can be cast in a more 
symmetric form with respect to the indices m, n, which 
we present in the main text in Eg. pip . 

We limit ourselves to superconductors with inversion 
symmetry. Then the singlet and triplet order parameters 
transform under inversion as follows, 

7'A(R, p) = A(-R, -p) - A(R, -p) = A(R, p) , 
7'A(R, p) = A(-R, -p) = A(R, -p) = - A(R, p) , 

where we assumed that the order parameter is an even 
function of the spatial coordinates R. This assumption is 
justified by the analysis of the behavior of the off-diagonal 
functions /m(p)- The expansion of the anomalous propa- 
gators in the Landau level basis contains all components 
(R|m), however, even and odd coefficients have different 
parity under inversion p ^ — p, 

/^,(-p) = E P™,„(|p|)A„(p) , 



'(p)2?m,„(|p|)A„(p) 



As a result, it is easy to show that for both singlet and 
triplet order parameters, the even and odd coefficients 
A„ are decoupled since no mixed term survives averaging 
over the Fermi surface. 



^X,(P)/-*(P). 



(B14) 



Note also that, in zero field for superconductors with 
basis functions (3^(p))fs = 0, the off-diagonal impurity 
self-energy vanishes since (/(p))fs = 0. Under magnetic 
field, however, the direction p is inequivalent to the per- 
pendicular to it direction, px, and /(p) is not simply 
proportional to 3^(p). Hence for the field in the plane 
and d-wavc gap (/(p))fs ^ 0, and there is a contribution 
to the off-diagonal self-energies from impurities. This in- 
tegral still vanishes for p-wave order parameters, since 
in magnetic field the directions p and — p may remain 
equivalent and so the symmetry /(— p) = — /(p) is still 
valid. 
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